source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/fdx1/bmsc/'

sobj <- read_rds('mission/fdx1/bmsc/bmsc.osteo.wtko.rds')

sobj <- sobj |> quick_process_seurat(skip_norm = T)

sobj |> write_rds('mission/fdx1/bmsc/bmsc.osteo.wtko.rds')

sobj |> DimPlot(group.by = 'fine_type', split.by = 'orig.ident')

# find all marker -----------
osteo_remarker <- sobj |>
  filter(orig.ident == 'WT') |>
  FindAllMarkers(logfc.threshold = 1, only.pos = T) |>
  as_tibble()

osteo_remarker |>
  filter(p_val_adj < .05, !str_detect(gene, '^Gm\\d+|\\d+Rik$')) |>
  slice_head(n = 100, by = cluster) |>
  write_source_csv('osteo_all_marker')

sobj <- sobj |>
  FindClusters(algorithm = 4, resolution = .5, random.seed = 1)

sobj |> DimPlot(label = T, label.box = T, label.size = 2, repel = T,
                cols = 'Paired') +
  theme_jpub(theme_classic)

publish_pdf('bmsc.osteo.leiden.umap.pdf', width = 60)

# determine cell type ----------
osteo_markers <- list(
  MSCs = c("Lepr", "Pdgfra", "Nes"),
  Cycling = c('Mki67','Top2a'),
  Pre_osteoblast = c("Runx2","Sp7", "Col1a1", "Alpl"),
  Osteoblast = c("Bglap", "Spp1", "Ibsp"),
  Chondrocyte = c('Col10a1','Acan'),
  Endothelial = c('Cdh5','Pecam1'),
  Osteoclast = c('Acp5','Ctsk')
)

sobj |>
  filter(orig.ident == 'WT') |>
  DotPlot(osteo_markers, cols = 'RdYlBu', cluster.idents = T, dot.scale = 4) +
  labs(x = 'Gene', y = 'Cluster') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('bmsc_osteo_subcluster_marker_dotplot',
                    width = 140)

sobj |>
  DimPlot(cols = 'Paired', label = T, label.box = T)

sobj <- sobj |>
  mutate(finetype = case_when(seurat_clusters %in% c(4,5) ~ 'Pdgfra+ MSC',
                              seurat_clusters == 6 ~ 'Bglap+ Osteoblast',
                              seurat_clusters == 9 ~ 'Runx2+ Pre_osteoblast',
                              seurat_clusters == 2 ~ 'Col10a1+ Chondrocyte',
                              seurat_clusters == 8 ~ 'Acp5+ Osteoclast',
                              .default = 'Cdh5+ Endothelial-like'))

# customer request ------------
sobj <- sobj |>
  mutate(finetype = case_when(seurat_clusters %in% c(4,5) ~ 'BMSC-1',
                              seurat_clusters %in% c(2,6) ~ 'BMSC-2',
                              seurat_clusters == 8 ~ 'Osteoclast',
                              .default = 'Pre-osteoblast'))

sobj |>
  DimPlot(group.by = 'finetype', split.by = 'orig.ident') +
  ggtitle('Prrx1+ Osteo-lineage') +
  theme_jpub(theme_classic)

publish_pdf('Prrx1_celltype_umap.pdf', width = 110)

sobj |>
  ggplot(aes(orig.ident, fill = finetype)) +
  geom_bar(position = 'fill') +
  labs(x = 'Group', y = 'Proportion', fill = 'Cell type') +
  theme_jpub()

publish_pdf('osteo.subtype.frac.pdf')

sobj |>
  DotPlot(osteo_markers, group.by = 'finetype', cols = 'RdYlBu', dot.scale = 5) +
  labs(x = 'Gene', y = 'Cell type') +
  theme_jpub() +
  RotatedAxis()

publish_source_plot('bmsc_osteo_annotated_marker_dotplot', width = 140)

# slingshot ----------
library(slingshot)

osteo_umap <- Embeddings(sobj, reduction = 'umap')

osteo_ss <- osteo_umap |>
  slingshot(sobj$fine_type)

osteo_pt <- osteo_ss |> slingPseudotime()

osteo_pt |> head()

sobj |>
  AddMetaData(osteo_pt) |>
  FeaturePlot('Lineage1') +
  ggtitle('Pseudotime: Lineage1')

sds <- as.SlingshotDataSet(osteo_ss)

sds |> glimpse()

sds@curves$Lineage1$s |>
  head()

sobj |>
  AddMetaData(osteo_pt) |>
  FeaturePlot('Lineage1', cols = c('lightgrey','red')) +
  ggtitle('Pseudotime: Lineage1 in Prrx1+ cells') +
  geom_path(data = sds@curves$Lineage1$s, aes(x = umap_1, y = umap_2),
            arrow = arrow(type = 'closed', length = unit(.1, 'inches')))

sobj |>
  AddMetaData(osteo_pt) |>
  FeaturePlot('Lineage2', cols = c('lightgrey','red')) +
  ggtitle('Pseudotime: Lineage2 in Prrx1+ cells') +
  geom_path(data = sds@curves$Lineage2$s, aes(x = umap_1, y = umap_2),
            arrow = arrow(type = 'closed', length = unit(.1, 'inches')))

## DDRTree ----------
shelf(DDRTree2)

data_mtx <- sobj |>
  GetAssayData()

data_mtx |> nose()

var_gene <- sobj |> VariableFeatures(nfeatures = 2000)

var_gene |> glimpse()

var_mtx <- data_mtx[var_gene,]

var_mtx |> glimpse()

FM <- var_mtx |>
  t() |>
  scale() |>
  t()

FM |> glimpse()

ddrt_center <- function(ncells, ncells_limit=100) {
  round(2 * ncells_limit * log(ncells)/(log(ncells) + log(ncells_limit)))
}

# 20s for 2000 genes x 5336 cells  
# 12s for 1000 genes x 5336 cells in openMP version!
system.time(ddrt <- FM |>
  DDRTree(ncenter = ddrt_center(ncol(FM)), verbose = T)
)

ddrt_lineage <- ddrt$Y |>
  t() |>
  set_colnames(c('ddrt1','ddrt2')) |>
  as_tibble()

ddrt_embed <- ddrt$Z |>
  t() |>
  set_colnames(c('ddrt1','ddrt2')) |>
  as_tibble() |>
  mutate(.cell = colnames(sobj))

sobj <- sobj |>
  left_join(ddrt_embed)

sobj |>
  ggplot(aes(ddrt1, ddrt2, color = fine_type)) +
  geom_point(data = ddrt_lineage, color = 'black') +
  geom_point(size = AutoPointSize(sobj@meta.data)) +
  theme_classic() +
  labs(x = 'Component 1', y = 'Component 2', color = 'Cell type') +
  guides(
    color = guide_legend(override.aes = list(size = 5))
  )

sobj |>
  ggplot(aes(ddrt1, ddrt2, color = finetype)) +
  geom_point(data = ddrt_lineage, color = 'black', size = .5) +
  geom_point(size = AutoPointSize(sobj@meta.data)) +
  labs(x = 'Component 1', y = 'Component 2', color = 'Cell type') +
  guides(
    color = guide_legend(override.aes = list(size = 5))
  ) +
  facet_wrap(~orig.ident) +
  theme_jpub(theme_classic)

publish_pdf('Prrx1_pseudotime_component12.pdf', width = 110)

## new marker set --------------
bmsc_marker <- c('Lepr','Adipoq','Lpl','Alpl','Col1a1','Col2a1','Acan','Col3a1','Postn','Plin1')

sobj |>
  DotPlot(bmsc_marker, cols = 'RdYlBu', cluster.idents = T,
          group.by = 'fine_type') +
  RotatedAxis()

sobj |>
  DotPlot(osteo_markers, cols = 'RdYlBu', cluster.idents = T,
          group.by = 'fine_type') +
  RotatedAxis()

# Id1/2 Slc31a1 ----------
sobj |>
  FeaturePlot(c('Slc31a1'), order = T, split.by = 'orig.ident') &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('Slc31a1_WT_KO_umap.pdf', width = 100)

sobj |>
  FeaturePlot(c('Slc31a1'), order = T,
              split.by = 'orig.ident', min.cutoff = 'q50') &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('Slc31a1_WT_KO_umap_cutoff.pdf', width = 100)

sobj |>
  FeaturePlot(c('Id1','Id2','Prrx1'), order = T,
              split.by = 'orig.ident') &
  theme_jpub(theme_classic) & NoLegend()

publish_pdf('osteo_Id1_2_prrx1_WT_KO_featureplot.pdf',
            width = 100, height = 150)

sobj |>
  filter(str_detect(finetype, 'MSC')) |>
  bill.violin(c('Id1','Id2','Slc31a1'), orig.ident) +
  labs(x = 'Group', fill = 'Group', y = 'Normalized expression',
       title = 'Prrx1+ Pdgfra+ Osteo-MSC') +
  stat_compare_means(comparisons = list(c('WT','KO')), label = 'p.signif',
                     vjust = .2)

osteo_kovwt_subdeg <- sobj |>
  FindMarkersAcrossVar(split.by = 'finetype', group.by = 'orig.ident',
                       ident.1 = 'KO')

# Col1a1 Ocn(Bglap) ---------
g1 <- sobj |>
  FeaturePlot(c('Col1a1'), order = T, min.cutoff = 'q10',
              split.by = 'orig.ident')

g2 <- sobj |>
  FeaturePlot(c('Bglap'), order = T,
              split.by = 'orig.ident')

g1 / g2

publish_pdf('Col1a1_Bglap_WT_KO_umap.pdf', width = 150, height = 150)

# DEG in BMSC1/2 -----------
osteo_sub_kovwt <- sobj |>
  FindMarkersAcrossVar(split.by = 'finetype', group.by = 'orig.ident',
                       ident.1 = 'KO')

osteo_sub_kovwt |>
  filter(gene %in% c('Id1','Id2','Slc31a1')) |>
  mutate(cluster = as.character(cluster)) |>
  ggplot(aes(fill = avg_log2FC, y = cluster, x = gene, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  scale_fill_gradient2(high = 'red', low = 'blue') +
  theme_jpub() +
  labs(y = 'Cell type', title = 'Prrx1+ cells: KO vs WT')

publish_source_plot('prrx1_id12_slc31a1_kovwt_logfc', width = 60)

osteo_sub_kovwt |>
  filter(str_starts(cluster, 'BMSC-1'), !str_ends(gene, 'Rik')) |>
  plot_pub_volc(group1 = 'KO', group2 = 'WT') +
  ggtitle('Differential expressed genes in BMSC-1')

publish_source_plot('bmsc1.ko.vs.wt.volcano', width = 60)

osteo_sub_kovwt |>
  filter(str_starts(cluster, 'BMSC-2'), !str_ends(gene, 'Rik')) |>
  plot_pub_volc(group1 = 'KO', group2 = 'WT') +
  ggtitle('Differential expressed genes in BMSC-2')

publish_source_plot('bmsc2.ko.vs.wt.volcano', width = 60)

## GSEA ------------
library(clusterProfiler)
bmsc1_gsego <- osteo_sub_kovwt |>
  filter(str_starts(cluster, 'BMSC-1'), p_val_adj < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db', keyType = 'SYMBOL', eps = 0,
        pvalueCutoff = .5)

bmsc1_gsekegg <- osteo_sub_kovwt |>
  filter(str_starts(cluster, 'BMSC'), p_val_adj < .05) |>
  batch_enrich_path(org = 'Mm', path = 'KEGG', method = 'GSEA')

bmsc1_gsekegg$`BMSC-1` |>
  plot_enrichment(padj_thres = 1) +
  theme_jpub() +
  labs(title = 'KEGG pathway enriched in BMSC-2: KO vs WT')

publish_source_plot('bmsc1.kovwt.kegg', width = 70)

set.seed(42)

bmsc1_gsekegg$`BMSC-2`@result |>
  mutate(p.adjust = pvalue) |>
  plot_enrichment(padj_thres = .7, force_regex = 'stem ce') +
  theme_jpub() +
  labs(title = 'KEGG pathway enriched in BMSC-2: KO vs WT')

publish_source_plot('bmsc2.kovwt.kegg', width = 70)

bmsc_kovwt <- sobj |>
  filter(str_detect(finetype, 'BMSC')) |>
  FindMarkers(group.by = 'orig.ident', ident.1 = 'KO') |>
  as_tibble(rownames = 'gene')

bmsc_gsekegg <- bmsc_kovwt |>
  mutate(cluster = 'BMSC') |>
  batch_enrich_path(org = 'Mm', path = 'KEGG', method = 'GSEA')

set.seed(42)

bmsc_gsekegg$BMSC |>
  plot_enrichment(padj_thres = .8, force_regex = 'stem ce') +
  theme_jpub() +
  labs(title = 'KEGG pathway enriched in BMSC: KO vs WT')

publish_source_plot('bmsc.kovwt.kegg', width = 70)

mm_msigc2 <- msigdbr::msigdbr(species = "Mus musculus", collection = "C2")

mm_msigc2$gs_subcollection |> table()

mm_msigc2 |>
  filter(str_detect(gs_name, 'PLURI')) |>
  distinct(gs_name)

mm_msigc2_t2g <- mm_msigc2 |>
  filter(str_detect(gs_subcollection, 'CP:WIKI')) |>
  dplyr::select(gs_name, gene_symbol)

bmsc1_gsec2 <- osteo_sub_kovwt |>
  filter(str_starts(cluster, 'BMSC-1'), p_val < .05) |>
  pull(avg_log2FC, name = 'gene') |>
  sort(decreasing = T) |>
  GSEA(eps = 0, TERM2GENE = mm_msigc2_t2g, pvalueCutoff = 1)

bmsc1_gsec2@result |>
  filter(str_detect(Description, 'STEM')) |>
  plot_enrichment(padj_thres = 1) +
  theme_jpub()

bglap2_gsekegg <- bglap2_deg$gene |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID', OrgDb = 'org.Mm.eg.db') |>
  left_join(bglap2_deg, join_by(SYMBOL == gene)) |>
  pull(avg_log2FC, name = 'ENTREZID') |>
  sort(decreasing = T) |>
  gseKEGG(organism = 'mmu', eps = 0, pvalueCutoff = .5)

bglap2_gsekegg@result |>
  plot_enrichment(metric = NES, padj_thres = .3) +
  labs(title = 'KEGG pathway enriched in Bglap2+ Pre-osteoblast: KO vs WT') +
  theme_jpub

bglap2_gsekegg@result |>
  write_source_csv('bglap2_gsea_kegg')

bglap2_gsekegg_res <-
  read_csv('mission/fdx1/bmsc/results/bglap2_gsea_kegg.csv')

## violin ------------
stem_plrp <- bglap2_gsekegg_res |>
  filter(str_detect(Description, 'pluripotency')) |>
  pull(core_enrichment) |>
  str_split_1('/') |>
  bitr(toType = 'SYMBOL', fromType = 'ENTREZID', OrgDb = 'org.Mm.eg.db') |>
  pull(SYMBOL)

select <- dplyr::select

sobj |>
  filter(str_starts(finetype, 'BMSC-1')) |>
  bill.violin(stem_plrp, group.by = orig.ident) +
  labs(x = '', fill = 'Group', y = 'Normalized expression',
       title = 'Pathway "Signaling pathways regulating pluripotency of stem cells"\nin BMSC-1') +
  theme_jpub(theme_classic)

publish_pdf('bmsc1.pluripotency.violin.pdf', width = 80)

sobj |>
  filter(str_starts(finetype, 'BMSC-2')) |>
  bill.violin(stem_plrp, group.by = orig.ident) +
  labs(x = '', fill = 'Group', y = 'Normalized expression',
       title = 'Pathway "Signaling pathways regulating pluripotency of stem cells"\nin BMSC-2') +
  theme_jpub(theme_classic)

publish_pdf('bmsc2.pluripotency.violin.pdf', width = 80)

sobj |>
  filter(str_starts(fine_type, 'Bglap2')) |>
  ScaleData(features = stem_plrp) |>
  DoHeatmap(stem_plrp, group.by = 'orig.ident')
